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A PROCEDURE FOR ESTIMATING STABILITY AND 
CONTROL PARAMETERS FROM FLIGHT TEST DATA BY USING 
MAXIMUM LIKELIHOOD METHODS EMPLOYING A 
REAL-TIME DIGITAL SYSTEM 

By Randall D. Grove, Roland L. Bowles, 
and Stanley C. May hew* 

Langley Research Center 

SUMMARY 

A maximum likelihood parameter estimation procedure and program were developed 
for the extraction of the stability and control derivatives of aircraft from flight test data. 
Nonlinear six-degree -of -freedom equations describe the aircraft dynamics and are used 
to derive the sensitivity equations for the method of quasilinearization. The maximum 
likelihood function with quasilinearization was used to derive the parameter change equa- 
tions, the covariance matrices for the parameters and measurement noise, and the per- 
formance index function. 

The maximum likelihood estimator was mechanized into an iterative estimation 
procedure utilizing a real-time digital computer and graphic display system. This pro- 
gram was developed for 8 measured state variables and 40 parameters. Test cases were 
conducted with pseudo or simulated data for validation of the estimation procedure and 
program. This program has been applied to a V/STOL tilt -wing aircraft, a military 
fighter airplane, and a light single -engine airplane. 

The appendixes describe in detail the particular nonlinear equations of motion, 
derivation of the sensitivity equations, addition of accelerations into the algorithm, opera- 
tional features of the real-time digital system, and test cases. 

INTRODUCTION 

The problem of estimating stability and control parameters of an aircraft from 
flight test data dates from the early days of flight. The results of early investigations 
were frequently limited, however, due primarily to insufficient estimation technology and 
restricted computational resources. Since 1960 there has been significant progress in 
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correcting these deficiencies; therefore, it seems appropriate to reevaluate, in a more 
general setting, the problem of estimating aerodynamic coefficients from flight data. 

Parameter estimation is the process of determining the parameters in a mathemat- 
ical model after having been supplied measured values for the variables of state and the 
input to the dynamic system. The accuracy of the resulting estimate is degraded by a 
combination of measurement, modeling, and numerical errors. Obtaining this estimate 
is a problem in inverse computation and the matter of existence and uniqueness of solu- 
tion must be resolved at least to some relative degree. The fact that in the defining of 
the discipline it is necessary to refer to modeling errors and the existence and xmique- 
ness of solution suggests that the expression "parameter estimation" is not adequate to 
describe the task, the task being a study undertaken by an analyst using the parameter 
estimation program. A meaningful consideration must include studies of the model 
definition in relation to the engineering application, the accuracy and precision to which 
the parameters are computed, and some indication of the uniqueness of solution. It is 
suggested that the expression "system identification" better represents the task being 
considered and more accurately implies the interrelated disciplines the study requires. 

A survey of the general problem of identifying a dynamic system from input -output mea- 
surements is given in reference 1. 

The objectives of this paper are twofold: 

(1) To present the development of an estimation procedure, based on maximum 

likelihood statistics, suitable for extracting stability and control param- 
eters of an aircraft from flight test data for realistic aircraft models 

(2) To develop a computer program and provide operational features of the estima- 

tion procedure when implemented on the Langley real-time simulation system 

The general approach adopted in this paper was based on a maximum likelihood output 
error method. The case where process noise is present is not considered. The assumed 
mathematical model for the aircraft was based on a standard six -degree -of -freedom rigid 
body description with linearized aerodynamic forces and moments. The interactive role 
of the analyst is discussed as well as various program options which are available. Also 
included in the paper is a demonstration of the performance of the estimation procedure 
and the computer program using pseudo flight data. The program developed has been 
successfully applied to the analysis of flight data for generically different aircraft 
(refs. 2 and 3). 
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SYMBOLS 


A sensitivity coefficient matrix (n X p') 

, 3 ^ 5^7 longitudinal, lateral, and vertical accelerations at center of gravity 

A, eg 1 jCg ^;Cg 

I’^ I’^Z I longitudinal, lateral, and vertical accelerations at instrument 

location 

B sensitivity coefficient matrix for N data points (nN X p') 

b wing span 

C^,Cm,Cn rolling-, pitching-, and yawing -moment coefficients 

Cl QjCn^o rolling -moment and yawing -moment coefficients at jS = 63^ = 6^. = 0 

Cm,o pitching-moment coefficient at «a = 6g = 0 

Cx>Cy,C 2 longitudinal-, lateral-, and normal-force coefficients 

^X o’^z o longitudinal -force and normal-force coefficients at = 6g = 0 

Cy,o lateral-force coefficient at ^ = 6^ = 6j. = 0 

c mean aerodynamic chord 

E mathematical expectation 

F vector function defined in equations of motion (n X 1 ) 

Fj components of F, where j = 1 , 2 , . . ., n 

G sensitivity equation matrix (n X n) 

g gravity 

H measurement noise vector for N data points (nN x 1 ) 

3 

L 



Ix^Y’^Z 

^XZ 

% 


Jn 

L 


^X’^Y’^Z 


m 

m' 

N 


n 

P’ 




Rl 

S 

T 

TxjTyjTz 

t 


ti 


u,v,w 


moment of inertia about X-, Y-, and Z-axis, respectively 

product of inertia 

wing tilt angle 

performance index function 

maximum likelihood function 

rolling, pitching, and yawing moments 
mass 

dimension of control deflection vector 

number of data points 

dimension of system 

number of parameters 

roll, pitch, and yaw angular velocities 

measurement noise covariance matrix (n x n) 

wing area 

flight test time period 

thrust along X-, Y-, and Z-axis, respectively 
time 

data point time, where i = 1, 2, . . ., N 
longitudinal, lateral, and vertical velocity components 
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V 

Vss 

Wj 

X 

a 




“i 


AQ! 

Aai 

6 

6-i 


" 1 ] 




P 



“i 





velocity 

slipstream velocity 
weighting matrix (nN x nN) 
state vector (n x 1) 

components of x, where k = 1, 2, . . n 

parameter vector (p' x 1) 

angles of attack and sideslip 

components of a, where i = 1, 2, . . p’ 

parameter change vector (p' X 1) 

components of Aa, where i = 1, 2, . . p’ 

control deflection vector (m* x 1) 

aileron, elevator, and rudder control deflections 

Kronecker delta 

measurement noise vector (n x 1) 

components of rf 

atmospheric density 

correlation coefficient of aj and cxj 

standard deviation of o!j 

roll, pitch, and yaw angles 
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Matrix exponents: 


T indicates transpose matrix operation 

-1 indicates inverse matrix operation 

The superscript M denotes measured value. 

The superscript o indicates nominal evaluation. 

Dot over a symbol denotes time derivative. 

Arrow over symbol indicates vector. 

BACKGROUND DISCUSSION 

Early methods for identifying aircraft stability and control parameters are typically 
characterized as equation error methods. They were basically least-squares estimators 
which minimize the equation error and are known to give biased estimates in the presence 
of measurement noise. Details of these various methods can be found in references 4 
to 7. Since the unknown parameters enter the equations of motion in a linear fashion, 
equation error methods are characterized computationally as single step processes and, 
therefore, are simple to deal with. 

More recent parameter identification methods are generally classified as output 
error techniques. These methods minimized the output error (measurement noise) 
between the measurements and the true states of the dynamic system and are often used 
to modify the initial estimates obtained by equation error methods. Unlike equation 
error methods, the identification problem using output error methods is nonlinear and 
this requires an iterative solution. Standard output error methods include the Newton- 
Raphson iteration method (ref. 8), modified Newton -Raphson method (ref. 9), quasilineari- 
zation method (refs. 10 to 14), and various forms of gradient -dependent methods (ref. 15). 
The quasilinearization and modified Newton -Raphson methods can be shown to be identi- 
cal. The Kalman filter, a sequential estimation method, can be shown for certain restric- 
tive cases to be equivalent to the two techniques just mentioned and, therefore, can be 
considered an output error method. Output error methods are known to produce unbiased 
parameter estimates under realistic conditions on the measurement noise. However, if 
a significant amount of process noise exists, that is, gusts and modeling errors, then the 
validity of estimates obtained by using output error methods is subject to serious 
question. 
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PROBLEM STATEMENT AND ASSUMPTIONS 


Assume that the equations of motion of an aircraft are in the form 

X = f(x ,c?,6’(t)) + ]l{t) 

with a measurement model defined as 
x^(t) = h(x(a,t)) + ^(t) 

6^(t) = 6(t) + f(t) 

The variables are defined as follows: 

X aircraft state vector 

6^(t) aircraft control vector 

a unknown parameter vector 

3T^ state measurement vector 

— M 

5 (t) control measurement vector 

Jl{t) process noise vector 

?f(t) measurement noise vector 

y(t) measurement noise on input 

h (x (a,t)) measurement output vector 

The essential problem is to estimate the parameter vector a when given the equa- 
tions of motion of the dynamic system and the measurements X (t) and 6 (t). 

The solution of the problem as posed in the preceding paragraph is probably not 
practical at the present time. Various approaches to the general problem have been 
attempted (for example, refs. 16 to 18). The results of these studies indicate that for 
dynamic systems as complicated as aircraft, a solution of the estimation problem 
requires substantial computational effort. In order to proceed with sC solution of the 
estimation problem which is computationally feasible and for which theoretical techni- 
ques have been developed, the following assumptions are made: 


8 


(1) Rigid body aircraft (six degrees of freedom) 

(2) Only measurement noise (fT(t) = r (t) = o) 

(3) The measurement noise v(h) is a sequence of independent Gaussian random 
variables with 

E = 0 

= Rl«ij 

where the matrix Rl is unknown. 

The foregoing assumptions imply the following conditions: 

(1) The aircraft maneuvers do not exceed a dynamic range consistent with lineariza- 
tion of the aerodynamic forces and moments. 

(2) Wind gusts, modeling errors, and inaccuracies of measuring physical movement 
of a control surface are considered sufficiently small to warrant the use of an output 
error method. 

(3) Alinement and location of rate gyros and angle of attack and sideslip instrumen- 
tation are of sufficient quality to preclude their inclusion in the measurement model. 
Anomalies introduced by accelerometer location are included in the measurement model. 

PARAMETER ESTIMATION PROCEDURE 

The parameter estimation procedure, using the method of maximum likelihood with 
quasilinearization, is diagramed in figure 1. The basic components are the (1) equations 
of motion, (2) sensitivity equations, (3) maximum likelihood estimation equations, (4) per- 
formance index, and (5) flight test data. These components are described next and the 
procedure in figure 1 is explained in a subsequent section. 

Equations of Motion 

The equations of motion are for a six -degree -of -freedom rigid body aircraft and 
are stated in detail in appendix A. The equations of motion are written in a general vec- 
tor form for simplicity in formulating the parameter estimation algorithm. 

The nonlinear equations of motion in vector form are 

X = F(x,a,6)= [ji. Eg, . . ., (1) 

where n is the dimension of the system (n ^ 12). 
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The state vector is defined as 

X = X (o!,t) = pij, Xg, . . .,Xi^ (2) 

where x (a,t) is the solution of the equations of motion and x{^,t) = x°(t) is the 
nominal solution. 

The parameter vector, the components of which are the system parameters, is 

T 

a = |q!^, a^, . . (3) 

where p' is the total number of system parameters and is the nominal parameter 

vector. These parameters are the aerodynamic coefficients (stability and control deriva- 
tives) and the initial conditions of the state. These parameters must be initialized for 
the first iteration of the algorithm. The coefficients are initialized by a prior estimate 
(wind-tunnel data or analysis) and the state is initialized from the flight test data. 

The input to the system is 

6 = 5 (t) = ®2’ ' * *’ ^m’J 

where S is the control deflection vector with dimension m'. The control vector 6 
is set equal to the measured control input to the aircraft 

6(ti)=5^(ti) (i= 1, 2, . . .,N) (5) 

with linear interpolation between points. 

Sensitivity Equations 

The sensitivity equations form a basis for the method of quasilinearization and are 
derived by formally differentiating the equations of motion with respect to the parameters. 
The sensitivity equations are integrated in parallel with the equations of motion to yield 
the sensitivity coefficients, which reflect the sensitivity of the nominal solution with 
respect to each parameter. The method of quasilinearization uses these coefficients to 
linearize the change in the solution of the nonlinear equations of motion due to a change in 
the system parameters. Reference 19 describes the use of these parameter sensitivity 
coefficients in dynamic systems. 

Sensitivity coefficients. - Used in the method of quasilinearization is a linear approx- 
imation expressing the change in the state vector as a linear function of the changes in the 
parameters. The expansion of the nominal solution about the nominal parameter vector, 
neglecting second and higher order terms, is 
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r 



where each partial derivative (sensitivity coefficient vector) is evaluated along the 
nominal solution. 

The equation in matrix form is 

3r(a°+Ao!,t) - jr(a°,t) = A(t) Aq! (7) 

where 



and 


Aa = [Aq!^, Aoig, • . 

This matrix equation expresses the change in the nominal solution as a linear function of 
the parameter changes and the sensitivity coefficients. This equation is used in the 
expansion of the maximum likelihood function about the nominal parameter vector. 

Derivation.- The sensitivity equations are derived from the equations of motion by 
taking the partial derivative of each equation with respect to each parameter. The sensi- 
tivity equations corresponding to the equations of motion in appendix A are derived in 
detail in appendix B. 

By assuming that the derivatives are continuous in a and t (ref. 20), 



that is, the order of differentiation can be interchanged. This result is used to derive 
the sensitivity equation 
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This equation in matrix form is 


_d 

dt 




+ 


8F 

3o!i 


where 


and 


G(t) = 



~9F 

9F 

9F~ 


8xj’ 

9 x 2 ’ • • 

” 9^n_ 


9Fi 

9Fi 


9Fj 

9xj 




9F2 



9F2 

9x]^ 

9X2 

• • • 

9Xn 

9F 

■^n 

9Fn 

. . • 

£n 

9X1 

9X2 


9Xn_ 

1 

CD 



T 

90!^^ ’ 

9Q!. ’ ' 

■’ 9Q!^ 



(i= 1, 2, . . .,p’) (9) 

(i = l, 2, . . .,p’) (10) 


This system of equations represents p' sets of n simultaneous first-order 
linear differential equations with time varying coefficients. The solutions of this sys- 
tem are the sensitivity coefficients where i = 1, 2, . . p' and 

k = 1, 2, . . n. These coefficients are the elements of the matrix A(t) used in the 
linear approximation of the change in the nominal solution. 

The initial conditions of the sensitivity coefficients corresponding to the aerody- 
namic parameters are 


and the initial conditions of the sensitivity coefficients corresponding to the initial con- 
ditions of the state (parameters) are 
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90!j^ 


9xi, 


t=0 


9x^(0) 


t=0 


= «ik 


where 6^^^^ is the Kronecker delta. 


Maximum Likelihood Estimation 

The maximum likelihood method is used to estimate the stability and control param- 
eters of the aircraft from flight test data. This method has the asymptotic properties of 
unbiased and minimum variance estimates (ref. 18). Maximum likelihood estimation 
requires initial parameter values to start the algorithm, but it assumes that the covari- 
ance matrix of the measurement noise is unknown. 

Maximization of the likelihood function yields the parameter change equations and 
the covariance matrix for the parameters. These equations yield the changes in the 
nominal parameters to improve the fit between the measured and calculated variables of 
state. The covariance matrix gives the variances (or standard deviations) of the param- 
eters, or in an inverse sense the sensitivity of the parameters in the equations of motion. 
This matrix also indicates the dependency or correlation among the parameters. 

Maximization of the likelihood function also yields the covariance matrix for the 
measurement noise based on the current nominal solution. This matrix gives the vari- 
ances (or standard deviations) of the difference of the measured state and the nominal 
solution. The inverse of this matrix is used in the parameter change equations as a 
weighting matrix. The performance index function is derived by substituting the covari- 
ance matrix of the measurement noise into the likelihood function (ref. 21). The per- 
formance index function derived is the determinant of the covariance matrix of the mea- 
surement noise which is defined as the criterion for the maximum likelihood method. 


Measurement noise.- Let the measurement noise be 


?7(ti) = x^(ti) - X (a,tj) 




(i = 1, 2, . . ., N) (11) 


where N is the number of data points, x^(tj) is the measured data, x (a,tj) is the 
calculated solution for the true a values, and r)(t^ are independent Gaussian vari- 
ables with zero mean. The noise is assumed to have the statistical properties of 



( 12 ) 



= 


where E denotes the mathematical expectation and 
matrix of the measurement noise. 


Rl 


(13) 

is the unknown covariance 
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Maximum likelihood function.- The maximum likelihood estimation of Aa and 

is obtained by maximizing the likelihood function L^a°+Ao!,Ri) with respect to Aa 
and R^, respectively. Although the term Aa is not explicit in the function L(a,Rj), 
its use is evident in the maximization procedure. 

The likelih ood function for the state is (ref. 18) 

N 


i=l 1 


(14) 


where jR^j denotes the determinant, R£^ denotes the inverse of the symmetric 

II _j|2 _*>p 

covariance matrix Rj, and l|iqL_l = x Rj^ x. 

Parameter change estim ation .- The maximum likelihood function for the nominal 
solution is 


N 

l(5e“,ri) = - 1 y ;ryti)Ryi7(ti)- f ^ |ri| 

i=l 


(15) 


where 


^(ti) = x^(ti) - x°(ti) 


and Rj^ is estimated from the nominal solution. The likelihood function is expressed 
for a change in the nominal solution by using equation (7) as 


x^(ti) - X (a°+AQ!,ti) = x^(ti) - x°(ti) - A(ti)Ao! 

= %) - A(ti)^ 


(16) 


Hence, 


L (a°+A^, Rl) = - I ^ { [%) - A(ti) Aa]"^ R-^[h (ti) - A(ti) Aa] | | In | Ri| (17) 

To maximize this equation, the partial derivatives of L(a°+Aa,R]_) with respect to each 
parameter change Aa^ are set to zero; that is. 


8L 

9Aa 


N ^ 

(a°+Aa,R^) = - - -2A'^(ti)Rj^rf(ti) + 2A’^(ti)Rj^A(t^)Aa 


= 0 


(18) 
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This operation implies that 


N 

i=l 


Vi 




( 19 ) 


In matrix form the equations are 


where 



Aa = b"^WjH 



A(tl) 


^(*l) 


Rl 0 

B = 

A(*2) 

H = 


II 

^1 


A(tN') 

(nN,p') 


(nN,l) 

1 

o 

1 


-1 


(nN,hN) 


( 20 ) 


For p' > n, there exists more unknowns than equations; the system is overdeter- 
mined by calculating the nominal solution and sensitivity coefficients at N data points. 
The solution of Ao; is given by 

Aa = (21) 

where the new nominal parameter vector is incremented by Aa. A necessary condition 

T 

for this to be the estimate to maximize the likelihood function is that the matrix B WfB 
must be positive definite. 


Maximum likelihood methods (ref. 18) give asymptotically unbiased estimates and 


( 22 ) 


the inverse of 
parameters: 

e[a^ = 0 

for 

e[h] = 0 

and 

e[aq! Aoi'^ = 

for 

eIjih"^ = Wj 


(23) 
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The error covariance matrix for the estimated parameters is (ref. 22) 




Oi^ 








CL€ 


^0! tQ'i‘^«p*^«l ^Q!p»Q!2^“p'^“2 




^“2“p' “2 “p’ 


“p’ 


( 24 ) 


2 

where a is the variance of the estimate, a is the standard deviation, and p 

“i “i“j 

is the correlation coefficient for and (Xj. 

Estimation of measurement noise statistics.- The likelihood function of equa- 
tion (17) for the change in the nominal solution is 


N ^ rp _ 

L(a°+Aa, Rl) = - I 2 [ \^{k) - A(ti)Aa - A(ti)Aa 

i=l ^ 


- ^ In I Ri 


(25) 


where Rj denotes the unknown covariance matrix. The estimated value of the covari- 
ance matrix Rl is denoted by Rj(N). 

In a similar maximization procedure as for Aa, the likelihood function is maxi- 
mized with respect to R^. The derivative with respect to a matrix is defined as follows: 


3Rl = ("1 ' f iil) 

The derivative of l(^+Aq!,Rj) with respect to Rj is set to zero. But (ref. 21) 

■^(«°+^,Rl) = 0 (26) 

is equivalent to 

-^(a°+^,Rl) = 0 (27) 

rv*r* “ 1 
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and, hence, 



( 28 ) 


or 

N rj, 

R’l(N) = - A(ti)AaJ - A(ti)AaJ (29) 

where Rj(N) is the predicted estimate of the covariance matrix for N data points due 
to the change in the nominal solution. 

The equation for calculating the estimate of the covariance matrix used in this 
algorithm is 

R°(N) = Estimate of R^ 


N 


i=l 

The schemes for calculating the estimate of Rj^ are discussed later. 
The matrix R^(N) is written as 

.2 


R°(N) 


CT 


^1 


^ 2^1 


riiV2 

u2 

^2 




^7iT?n 


CT, 


’?2^n 


V, 


n 


(30) 


(31) 


is the covariance of and rjy 

Per formance index . - The performance index or index function evaluation gives a 
measure of performance for the iterative estimation procedure. Selection of the index 
function for the maximum likelihood estimator is an important condition as to whether the 
estimation procedure converges to the true parameter values. (See ref. 21.) 


where is the variance of 




and a. 


Vi 
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The index function is derived from the likelihood function (ref. 21). Substituting 
the covariance matrix estimate R°(N) (eq. (30)) into the likelihood function 
(eq. (15)) gives 


L(a°,R°(N)) = - |(nN) - In 



(32) 


Thus maximization of the likelihood function is equivalent to the minimization of 



which is defined as the index function for the maximum likelihood estimation procedure. 
The minimization of R'j^(N) with respect to Aa yields equivalent parameter change 
equations as in equation (19). 


Flight Test Data 

The flight test data are composed of the onboard instrument measurements of the 
aircraft behavior and are assumed to be the output of the aircraft mathematical model 
superimposed with instrument noise. These data contain many individual aircraft maneu- 
vers stored on one magnetic tape, with each maneuver easily accessible to the central 
memory of the computer. These data are used for comparison with the mathematical 
model output and for initialization of and control input to the equations of motion. The 
measurements, x’^(ti) and 5 (tj) for i = 1, 2, . . ., N, are known for all states and 
control deflections corresponding to the equations of motion. 

Steps in Procedure 

The steps in the maximum likelihood estimation procedure, corresponding to fig- 
ure 1, are as follows: 

(1) Initialize the system parameters, where denotes the nominal or current 

values of the parameters. 

(2) Integrate the equations of motion and the sensitivity equations to obtain the 
nominal solution and the sensitivity coefficient matrix, respectively. 

(3) Form the comparisons of the flight test data and nominal solution for each data 
point time tj, where i = 1, 2, . . ., N and tj^ = 0 and tj;f = T. 
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(4) Form the maximum likelihood estimation equations from the comparisons in 
step (3) and the sensitivity coefficient matrix in step (2) (dash lines indicate accumula- 
tion of information over the flight test time period T). 

(5) Calculate the performance index 

(6) Calculate the parameter changes Aa and the statistical information matrix 
RJ(N). 

(7) Update the nominal parameter values in step (1) to start the next iteration pro- 
cedure and repeat steps until convergence. 

Each iteration of the procedure extends over the flight test time period T and 
results in the update of the parameters. Evaluations within the period are at specified 
data point times tj, where the intervals tj^j^ - tj are integer multiples of the integra- 
tion step size. The integration step size is made compatible with the flight data intervals 
and the problem dynamics. 

Acceleration measurements and equations were added later to improve the estima- 
tion procedure and the addition is presented in appendix C. 

COMPUTATIONAL CONSIDERATIONS 

The flight tests do not always necessitate the use of all the state variables and 
parameters in the estimation algorithm (n state variables and p’ parameters) for 
specific cases. These cases involve only a specific portion of the program, as with an 
excitation of only the longitudinal motion of the aircraft. The total number of differen- 
tial equations in the algorithm for evaluation and integration is n + np’; in appendixes A 
and B, n = 8 and p’ =40. 

The parameter estimation algorithm was programed with a variable dimensioning 
capability with respect to the number of state variables and parameters necessary for 
any specific case. The analyst through the operational control features (appendix D) 
could select any subsets of the state variables and parameters to be active in the param- 
eter estimation algorithm. Computer program parameters were activated for each state 
variable and each parameter desired. This operation generated two sequences of num- 
bers specifying the state variables and parameters which were active in the algorithm 
and neglected the remaining ones. 

This variable dimensioning of the algorithm allowed flexibility in the parameter 
identification study in that the analyst could alter the program easily for each specific 
aircraft maneuver or computer run. In addition, the number of integrations was reduced. 
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In the calculation of Aa the matrix equations (eq. (20)) were 

^ = b’^w^h 

If the matrix B and the vector H were calculated for N data points, this would 
result in the storage of an nN x p' matrix and an nN vector. The storage would 
increase as the number of data points N increased. To eliminate the storage being 
dependent on N, equation (19), which is restated, was used: 


N I'l 

i=l 


N 

I 

i=l 


where 


N 

b'^W^B = ^ A'^(ti)R'^A(ti) 

i=l 


and 


N 

b'^WjH=^ A'^(tj)R‘V(ti) 
i=l 


(34) 


(35) 


In these equations the matrix products are formed as a function of time and the dimen- 
sions of the matrix products were not a function of the number of data points N. In 
fact, the dimensions depended only on p' and n, the number of parameters and state 
variables. 

In the calculation of the measurement noise covariance matrix, three computational 
schemes can be used. The matrix can be updated (1) on the same iteration as Aa, 

(2) one iteration behind and (3) with a two step procedure. 

Scheme (1) uses equation (29): 

N rj, 

Ri(N) = ^ ^[%) - A(ti)AQ!][%) - A(ti)^] 
i=l 

In this equation the matrix A(ti) and the vector ^(t^) must be stored for each incre- 
ment of time until Ao! is calculated. 

Scheme (2) is easy to incorporate into the program by using equation (30) : 

Rj(N) ^ Estimate of Rj^ 


N 


i=l 
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In this equation the covariance matrix is calculated from the nominal solution and not 
with the predicted change in the nominal solution. The matrix is one iteration behind in 
the algorithm but the effect is negligible when the change in the solution is small, that is, 
when convergence is achieved. 


Scheme (3) uses a two step process for each parameter update. The first step is 
to calculate Aa in the usual manner. The second step is to update the parameters and 
integrate only the equations of motion. Thus, the covariance matrix is 


N 


i=l 




X (a°+Ao! 


j^i) 


^^(ti) 


"(“°+A0!,ti)J 


(36) 


Scheme (3) is similar to scheme (1) in that scheme (1) uses the predicted change, whereas 
scheme (3) uses the calculated change in the nominal solution. Scheme (3) approaches 
scheme (2) when Aa approaches zero. 

Scheme (2) was used in the parameter estimation program with the option of using 
scheme (3). In test cases using schemes (2) and (3), the indication was that the differ- 
ence is not significant. 


TEST PROCEDURE 

The testing procedure used pseudo flight data in checking the maximum likelihood 
estimation algorithm. The data were generated by integrating the equations of motion 
and then adding measurement noise, all states assumed being measured. The measure- 
ment noise was sequences of pseudorandom numbers (random within the capability of a 
digital computer) having the normal (Gaussian) distribution with zero mean and known 
standard deviation. These data were assumed to be the flight test data or measured data 
for specified parameter values. The parameter values were then offset to become the 
nominal parameter values for the parameter estimation algorithm. 

Test cases using the pseudo flight data were conducted for the longitudinal motion 
of the aircraft and are presented in appendix E. The maximum likelihood algorithm com- 
puted the standard deviation of the measurement noise and the parameter values and 
their standard deviations; no statistical information was assumed concerning the noise. 
Results obtained from the test cases indicated that the calculated parameter values and 
standard deviations of the noise were converging to the true values. 

CONCLUDING REMARKS 

A maximum likelihood parameter estimation procedure and program have been 
developed and validated for the extraction of the stability and control derivatives of 
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aircraft from flight test data. A nonlinear six -degree -of -freedom aircraft mathematical 
model was used in the derivation of the sensitivity equations. Instrument measurement 
noise was accounted for by the maximum likelihood estimator. The program was devel- 
oped for 8 measured state variables and 40 parameters, from which subsets could be 
selected for program operation. Real-time digital simulation and graphic display pro- 
vided the analyst with interactive control and display capabilities during the study. The 
program has been applied to a V/STOL tilt -wing aircraft, a military fighter airplane, and 
a light single -engine airplane. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., April 6, 1972, 
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APPENDIX A 


EQUATIONS OF MOTION 


The equations of motion are stated for reference and in particular for the deriva- 
tion of the sensitivity equations in appendix B. The following equations of motion are for 
a V/STOL tilt -wing aircraft: 


u = -qw + rv - g sin e + 1 ^ + % H) + 


V = -ru + 


pw + g cos 0 sin 0 + I ^ V^S (cy^o + ||: + 


pb 
P 2V 


(Al) 


+ ^Yr ^ 


(A2) 


w 


= -pv + qu + g cos 0 cos (^ + I ^ V^S (^2,^0 


Olf 


“a 


+ Czq + CZ5g5e) + 


m 


(A3) 


p= qrl 


. |^(pq . r) . I ^ v2sh L . q . § 


+ qg^SrJ + I f-(yss + vysh(ci^^ cos i^ - sin iw\6a + J 


M- 


X 


^X 


^ ^ 2 # v2ss(c„,o * . C„ 


2V 


+ ^m,~ ^m;; 5 c I + 


q 2V ly 




(A4) 


(A5) 


r = pqf^^^ + ^(i, - qr) -H I V^Sbk^o + • § h- 


iz 


pb_ 
p 2V 


+ Cnr I7 + Cns/rjl - | l|(VsS + cos ^ (A6) 
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and 


e = q cos 0 - r sin (f> 

0 = p + (q sin (f) + r cos 0)tan 9 




q sin (j> + r cos 0 
cos 9 


V = \u^ + 


Q!a = tan 


-1 


w 

u 



|3 


= sin 


-1 


v_ 

V 



> 



(A7) 

(A8) 


(A9) 


Trim conditions of the aircraft are added by substituting, respectively, for /3, 

5a> 5e, and 6r in equations (Al) to (A6) the following terms: 

«a - «a,t 

®a - 
" ®e,t 

6r - 6j,^t 

where the subscript t denotes the trim conditions (values) of the aircraft. 

The equations of motion are altered in two ways: (1) solve for v (eq. (A2)) 

explicitly, and (2) decouple the p (eq. (A4)) and r (eq. (A6)) equations. The equations 
of motion are then written in simplified notation for use in the derivation of the sensitivity 
equations. 

The equations of motion are 

X = F(x,a,6,V,aa,aa»/3,/3) 

= [Fi, F2, . . Fg]"^ (AlO) 
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The state vector is 

X = |Xj, Xg, . . 

= [u,v,w,p,q,r,0,<^'^ 

The parameter vector is 

T 

a = jo-p a^, . . ^4^ 

= ^(0), Cx^o’ ■ • •’ 

The control deflection vector is 

® “ pa>®e>^^ 

Equations (Al) through (A8) are written as 
u= Fi(x,a,6,V,a!a) 

Tx 

= -qw + rv - g sin 0 + ajV2(Cxi + Cx2) + ^ 

Tx 

= -qw + rv - g sin 0 + aiV2^Cx,o + + 15T 

V = F2(x,a,6,V,|3) 

2 / \ 

-ru + pw + g cos 0 sin 0 + aj^V f Cyi + Cy 2) + 

^ 

T 

-ru + pw + g cos 0 sin 0 + aiV^^Cy^o ^Y^^ + + ^3^ (^YpP + ^Yr^) liT 

1 - agCy^ 




2 \ 

= -pv + qu + g cos 0 cos 0 + aj^V ^C 2 i + ^ZZ) * liT 
= -pv + qu + g cos 0 cos 0 + a^V^(^C2^o + 



(All) 

(A12) 

(A13) 

(A14) 

(A15) 

(A16) 
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APPENDIX A — Continued 


P = F4(x,a,6,V,/3,4) 

= agF4(x,a,6,V,^,^)+ bjFg(x,a,6,V,^,^) 

= agjbgqr + IxzPq + + C^2) + 

+ bjjbgpq - Ixz^l^ + ^4V^(Cnl + Cj^2) " ^4Vs^(Cns) + 

= agjbgqr + IxzPq + ^4V^(Ci^o + + C^^rj 

+ + Mxj + bijbgpq - Ixzqr + ^4y^{Cj,^o + 

+ a5V(c„^|3 + CnpP + Cn^r) - a4Vg2 (c;,g^6a) + 

q = F5(x,a,5,V,Q!a,aa) 

= a^pr + b4(r2 - p2) + a3v2(Cmi + C^2) + ^ 

= a^pr + b4(r2 - p2) + agV2^Cm,o + ^rua^^a. + 

+ agV^Cm^^aa + Cm^q^ + — 

r = Fg(x,a,6,V,^,/3) 

= h^F'^(x,a,d,V,p/) + bgF^(x ,a,6,V,i3,4) 

9 = Frjix) 

= q cos (p - r sin cp 
<p = Fg(x) 

= p + (q sin (p + r cos 0)tan 9 


(A17) 


(A18) 


(A19) 


(A20) 


(A21) 
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where 


a-^^c! 
1 " 2 m ® 

aj = aig) 

^4 = |pSb 

II 

Irt 

o _ ^z ■ ^X 
7 T 

a - 1 P a;; 

bj - ^xz 

Vz “ ^xz 

b2 = Iy - Iz 

h -^XZ 
°4 T 

Iy 

bg= ?X. 

^X^Z ^ 

and 


^S = Vgg + V 


^X1 = ^X,o + ^Xq, “a 



'^Yl = Cy^o + + CYg^6r 

Czi = Cz^o + C2^^o-a+Czg^ae 


^1,0+ + CiQ^^r 




^nl C„^o + + Cng 6j. 


" %a ^^6a % 


a3=aif^ 


ag = 


II - T^ 
X^Z ^XZ 


^9 = "8 J 


*^3 - Ix - h 


^X2 = % H 


C _ t) 

Cy2-2v 


(%pP + Cy^r^ 


Cz2 = Cz 

^q 2V 


Cz2-^(Cz^^ + qpP + q^rj 


cos i^ - Cng^ sin i^ 

^^3 - C^g^Sa 

■'■ ^“io!a“a + Cmg^Sg 

»H 

'“1^ 

II 

CM 

3 

u 




Cn3 - ^ng 6a 


APPENDIX B 


nF.T ATT. fi OF DERIVATION OF SENSITIVITY EQUATIONS 


The sensitivity equations for the method of quasilinearization are derived in detail 
for the equations of motion in appendix A. 

The sensitivity equations are 

d / 9 x\ ^ / 95 T\ 
dtUo!i) [dai 


_ I Y 3 F 3 F Y 

“ \Z 9Xk\9ai/ 9V Z 


^ 8 V ^ 8 F Y 


9 xk\aaiy 


3 F ^“a 


9Q-a 9^1 v^y 9^ 


8 F 

8|3 j 

^ 9 X 2 \ 

8 ^ 

9 X 2 * 



8 V 


k=l 


9 xk\ 9 o'i 


8xv 


3 F ^ y ay / 9 Xk\ 

a^ 9 V z^ 9 Xk\ 9 o;iy 


8 F ^^aj 

fdx^\ 

8 o!a 8 x 3 ' 

\daij 


8 F 

Ba\ 




8 F ^®a| 

^ 3 X 3 \ 

8 <ia 8 x 3 

\dai) 


8 F 8/3 
8|3 9X2 V“i> 


8 F 

80!^ 


where 

G'(t) = jgjk(t^ (h k = l, 2 , . . 8 ) 

The functions F 2 and F 3 do not contain o-g^ or / 3 - Thus, 

8 

/ nv \ flV' 8 q!o 

G' 


d / 3 x\ 
dtUofij " 


,,,.78x \ 8F ^“a Y 0-’ /"^Xk^ 8F 8g Y p’ 

^ \ 9 «i/ 9 d!g 8x3 Z 3 ky 8 Q!j^y 0|3 8X2^ 2k\^80!j 

k~X k— 1 


(Bl) 






? 80 !., 8 F 


8 F 


8Fr 


a 3 ^ 8 F ag 2 


Yq!^ 8 afg 8 x 3 9 Q!i 9 ^ 8 X 2 ^“i 




(i = 1, 2, . . 40) 


(B2) 
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where 


8X 

80!j 


~8u 8y 8w 9 p 8q 8r 80 80 ~ 

8o!i’ dai 8o!i’ Bai 8o!i’ aofi’ Soti’ Soti 


iT 


G(t) = |g.k(t)j (i, k= 1, 2, . . 8) 

-► T 

^(“i) =Pl(“i)’ ^2(“i)> • • •’ 

The following equations are used in the derivation of the sensitivity equations: 


(1-8) 

(1-6) 

(1,3,5) 

(5) 

(2,4,6) 

(2,4,6) 


8F3 av ^ 

8Fj 90 !^ 

9F. 

80 9 V 

1 8 Xk 

8V 9Xk 


ada 9xk 

8^ 8xj^ 

80 8V 9 xjj 


(4,6) (5) (4,6) 


8Fj 8V \ ^ 

8^ 8V 


da^ 8x3 8/3 9X2 



! ^Fj 9aa QFj ag 

8(ig^ 9x3 ^ 3 k ax2 

(1-6) (5) (4,6) 

~ 8 Fi 9F-j ada 9 F 3 8F4 qo SFg 

fry A = , J 4 . — J — + A — c- 

i( 8o!j^ adg^ 8x3 9 o!j^ Qp 3x2 ^*^i 


(i, k= 1, 2, . . ., 8) (B3) 




(B4) 


The numbers in the parentheses above each term indicate the derivations in which they 
are used; this is in reference to the equations of motion used in the derivation of sensitiv- 
ity equations. 


(1) Sensitivity equations derived from u equation (eq. (A14)): 



(i= 1, 2, 


., 40) (B5) 
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where 

gll = aiu(2Cxi + Cx2) - aiV2Cx^^ 

gl2 = r+aiv(2Cxi + Cx2) 

§13 = ^l'^(2Cxi + ^x2) + ^iV^Cxq, 2 ^2 

** u + w 

§14 = 0 


§15 


= -w + aoVC 


2''"-Xq 


§16 = ^ 


gl 7 = -g cos d 


§18 


= 0 


and 


3Fi 


= a,V^ 


^^X,o ^ 


= a. V^q;^ 

aCx ^ ^ 


8Fi 


Q!e 


dC 


Xq 


= a„Vq 


(2) Sensitivity equations derived from v equation (eq. (A15)): 

8 

/ HXi _\ 9 

+ ^ (i = 1, 2, . . ., 40) 


d / 9v\ _ 


dt V9o!i) ^2ky 9Q!^ j ' 9Q!i 


(B6) 


where 


-r + aj^u 


§21 = 


2Cyi + Cy2 - Cyov(u 2 + w^) 


- 1/2 


1 - agCy^ 


aiv 


^22 


2Cyi + Cy2 - Cy V (u^ + w^) + ajV2Cy„(u2 + w^) 


p + a^w 


^23 = 


1 - ^3Cy^ 

2Cyi + Cy2 - Cy^v (u2 + w2) 


^ ■ ^3^Y; 
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and 


S24 = 

^25 
^26 "" 


^27 " 


^28 “ 


w + agVCYp 


1 

“ ^3^Y^ 

0 


-u 

+ agVCy^ 

1 

- ^3^Y^ 

-g 

sin 6 sin <p 

1 


g cos 6 cos <p 

1 

" ^3^Y^ 


8F2 


8F2 

9 Cy,o ^ 

■ ^ 3 ^Y^ 

9 Cy^ 

8F2 _ 

a3Vp 

8F2 

9 Cyp 1 

- agCy . 

9 Cy^ 


a^v2)3 

8 F 2 _ 

a3V 

^ ~ ^3^Y^ 

8 Cy . 

/3 

1 - agCy^ 

a 3 Vr 

8 F 2 

ajV^Sr 

1 - aoC-y • 

3 3t(3 

8 Cy 

Or 

1 - ^3^y^ 


(3) Sensitivity equations derived from w 


equation (eq. (A16)): 


dt 




+ 


8^ 

8Q!j 


(i= 1, 2, . 


40) 


(B7) 


where 


g 31 - q + a^u(2C2i + C22) - ~2 2 

“ u + w 

S32 = -P + ^iv(2C2i + C22) 

§33 = a^w(2Czi + Cz2) + 


^34 “ 


= -V 


ggg = u + agVCzq 
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^36 “ ° 

= -g sin e cos (f) 
^38 " ^ ^ 


and 


=ajV2 


9Cz,o 


= a. v2q! 


9Fo 

^=“2''q 


8F 


3 O 17-2« 

= aiV 0, 


9Cz. ■‘^1’' 


(4) Sensitivity equations derived from p equation (eq. (A17)): 
p= F4(x,a,6',V,^,/§) 

= agF4(x ,a,6 + bjFg(£,a,6 ,V,^,/3) 

Equation (B3) for becomes 

= s, I ^ J§- + ^ 9V ^ 

&4k av 9Xk 9/3 9Xk 8/3 9V Sxjj 9V 9xj5-, 


'dFa dF'a 


I w 6 9V 9/3 9/3 9V 9g 9V 

ly9xk 9V 9Xk 9/3 9Xk 8^ 8V 9Xk 8V 9Xky 


„ ^ „ I as ^ jA 

=‘6B4k ^ bjgek -S^jsgk 




= S4k ^ -^|g2k 


and equation (B4) for becomes 


/9F4 8F4 .n dF^ fdFa 9 Fr a 4 

F 4 K) = a„ -^ + -4 ^ + bJ^ + — ^ ^ — 2 
4V 1 ; 6^90,. gx2 Bcx^l l\^ 80 !i g^ 9X2 ®“i 


(B8) 


= ^6^4(“i) + '^1^6(“i) 


(B9) 
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Then, 




Air + ^rr ■:x Sok ) rsr— + 1 a.% 


sf'a 9F!. aA 9F, 


k=l 


^4k 


8/3 9v 


80 !i 8/3 8 v 90!i^ 


+ b. 


Y / „ ^ \ ( Qxk\ / ^ 

/j I ^ 6 k 8 v ^^j\&oii) \9o!j^ 0 | 8 v 90!jy 

k~l ' 


dFl 


Z Kb + =>6 ^ g2k + >=1 av Kr 
k— 1 \ 

+ ^6^4(“i) 

Kxk 


8 


S4k a^rF4(“i) 


(i = 1, 2, . . ., 40) 


k=l 


where 


§41 = ^4^ 

^Z2 ^^Z3 

§42 = ^4^^ 

2Ca + C^2 + ^^Z3"^ ■ 


„2^ / 2 2 )"^/^ 

+ 

V + w^j 

§43 = ^4^" 

2Ca + C^2 + ^^13 


S44 = 


2)-1/2 
+ w^j 


^/3 2V 




+ w‘ 


\-l /2 



S 45 = ^ 2 ^ ^XZP 


S46 = V + 


IT IT r\ 

S 47 - ^48 " ° 


(BIO) 
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2Cnl + Cn2 - 2Cn3 " Cn^v(u2 + _ Cn^ 


§62 = 




2Cni + Cn2 - 20^3 ^ - Cn v(u2 + _ n |b 


V 


i3 W 


a 4 v 2 Cn.(u 2 +w 2)“^/2 


+ a 4 V“Cn^ 
§63 = a 4 W 


2 Cni + Cn 2 - 2 Cn 3 ^ ~ Cn^v(u 2 + w 2 ) - Cn^ ^ 


^64 "" ^ 3 ^ + 
^65 " 

^66 "%Z^ + 

^67 " ^68 = ® 


9F4 

_ "Fg _ 

^^1,0 

^^n,o 

9 F 4 

9F6 

9C;. ■ 

^)3 

ac„.-5 

9 F 4 

. < 

9C; 

tr 

9Cnr ' ^ 

9 F 4 

9F^ 


®^”5a 

8Fj 

3F^ _ 

«Cn 6 ^ 

^^Z6a 




= -a^Vg 5 ^ sin 


9)3 av ® '■p 


9 F' 8 F^ 



II 

oa 

a 

0 

CD 

II 

9 F 4 

9 Fe 

9 C, 

ip 

9 Cnp 

0 F 4 

9 F' 




= a 4 v 26 r 


9 Fe ao 

iM- a r 

9^ 9v'^5% 
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(5) Sensitivity equations derived from q equation (eq. (A18)): 

V • 9 “a /9F5 9Fg aPg' 


±(^\=SL' 
dt Z'A ^“a ^ 


+ "T" + TT 


ao!i y \9Q!i 9o!a 9w ao!j 


where 


= 1 ®5kfe 


k=l 


^5(“i) 


(i = 1, 2, . . .,40) 


g'51 - agu(2C^l + Cjn2) - ^^2 +^2 " 5 

^52 ~ *^ 1 ^ 2 ) 

g'53 = a8^(2Cml + Cm2) + ^sV^C^Q-a -2“7^ 

g*54 = a,^r - 2b4P 

4 5 = 

§56 ^?P ^^^4^ 


S57 “ ^58 - 0 


and 


2 


9F. 


9Fp; 

ar-^ = agVd-. 
aCjtn - « 


Oil 


and 


aFc 


9Cmf 


aoVq 


aF 


ac 


^=agV26e 


m5e 


a^a 9w a 


(Bll) 
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(6) Sensitivity equations derived from r equation (eq. (A19)): 

d fdr '' 
dt 


■) = h 



Lk=i\ 


Y / M ^^’4 ^ V 9xk\ /®F4 ^ ^^4 dB ^^2^ 
Z^|^^4k+ dv^^^j\dai) \^9Q!i dvdaj 

vfr, 

-1\ ^ 8v^^^\^9o!iy l^Q^i 3j3 9v 8a^J 


v/' ^ ^^4 80 ^ ^^6 80 \f^^k 

= Z ( ^6k av ^2k + ^ 3 ^ g2k 


axk^ 


3/3 


k=l\ 

+ V4(“i) + VeK) 


8/3 9V 




I - ?6(“i) 

k=l ' ' 


(i = 1, 2, . . ., 40) (B12) 


where all the terms have been defined in the derivation of the sensitivity equations for 
p equation. 

(7) Sensitivity equations derived from 6 equation (eq. (A20)) : 



(i = 1, 2, . . ., 40) (B13) 


where 

^71 ^72 "" ®73 = §74 ^ 


g -75 = cos 0 
§76 " ^ 

§77 = 0 

§78 ^ ^ 

(8) Sensitivity equations derived from 0 equation (eq. (A21)): 



(i = 1, 2, . . ., 40) (B14) 
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where 

^81 = ^82 = ^83 = 0 

^84 = ^ 

ggg = sin (p tan 9 
ggg = cos p tan 6 

_ p 

^87 cos 6 
ggg = 6 tan 6 
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ADDITION OF ACCELERATIONS INTO ALGORITHM 


The acceleration measurements and equations were included in the parameter esti- 
mation algorithm to improve the extraction process. They are used with or can replace 
the linear velocities u, v, and w. The acceleration equations were transformed to the 
instrument location from the center of gravity. 


The equations for the center of gravity are 


^X,cg 


u + qw - rv + g sin 6 

^Y,cg 

_ 1 
g 

V + ru - pw - g cos 9 sin (j> 

_^Z,cg_ 


w + pv - qu - g cos 6 cos <p 


(Cl) 


The equations are transformed to the instrument location (ref. 23), that is, 


^,I 


^,cg 


^,I 

= 

^,cg 

4 

^Z,I_ 


_^Z,cg_ 



(pq + r) -(p2 + r2) (qr - p) 

(pr - q) (qr + f>) -(p2 + q2) 


X. 


(C2) 


where are the center -of -gravity offsets of the accelerometer measurements. 

The acceleration sensitivity equations were derived in terms of the sensitivity 
equations and coefficients stated in appendix B and need only to be evaluated and not 
integrated. 

(1) Sensitivity equations derived from a^ j equation: 



+ 



+ z.. 




(i= 1, 2, . . ., 40) 


(C3) 
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(2) Sensitivity equations derived from Oy j equation; 

■(■Si) - "(■^) - 2yaP)(Si)+ Kp * "4^, 


^2i = i| 

90!j g 


(u - 2y^r + 2a<l)('gT) * S 6 sin « (-gr) ‘ g =»= « ““ 


9v 

dai 


— Zf2 


80!i 




da^ 


(i = 1, 2 40) 


(C4) 


(3) Sensitivity equations derived from a^ j equation: 


9a 


2^1 = 1 


90!j^ g 


■“(Si) ♦ “(Si) + - "“g“)(S^+ (■" + - 2aa<l)(-^ 


+ (*aP + ya<l)('^') + g ain e cos •#>(^) + g cos 0 sin 


\9«i 


\daij 


+ 



+ y£ 




ayao!. 



(i = 1, 2, . . .,40) (C5) 


The maximum likelihood function was modified to include the accelerations in the 
algorithm. The likelihood function is 


where 


N 

L(a^R2) = - I 2 | ^ |R2 


(C6) 


V (ti) = 


77 (ti) 

Hk)_ 



a^^l(ti) - a^j(ti) 
■ ^Z,l(*i) 
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APPENDIX C - Concluded 


and R2 is measurement noise covariance matrix with accelerations included; that is, 
RgCN) ^ Estimate of R2 

i=l 

The maximization procedure is similar to the previous developments and similar esti- 
mation equations can be derived. 
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APPENDIX D 


OPERATIONAL FEATURES 

The computer program was written in FORTRAN IV language (75 000 octal loca- 
tions) and run on the Control Data series 6000 digital computer complex, a major applica- 
tion being real-time simulation (RTS) (ref. 24). Incorporated in the RTS system are the 
cathode ray tube (CRT) graphic display units. 

The computer program was mechanized into an iterative estimation procedure with 
manual interactive control through the utilization of the RTS system. The operational 
diagram of the RTS system is shown in figure 2, the main components being the computer 
complex, control console, and CRT. The remaining equipment is for output of informa- 
tion and monitoring the program. Figure 3(a) shows a photograph of the program control 
station and figure 3(b) shows a closeup of the control console. 

The maximum likelihood estimation program resides in the central memory of the 
computer. The analyst investigating the stability and control derivatives of the aircraft 
has direct control of the computer program through the control console. The control 
console has mode control switches for program operation, a data entry keyboard for 
inputing program parameters and logic controls, logical switches for program options, 
and indicator lights for program status. The digital decimal display was used to monitor 
continuously any selected parameter or variable in the program, particularly the perfor- 
mance index function. 

The CRT displayed the flight test maneuver at the start of each iteration. The 
response of the equations of motion was plotted simultaneously as it was computed in the 
digital program and was plotted with the flight test maneuver for direct comparison. 

This display permitted quick analysis of each flight test case on an iteration to iteration 
basis. Figure 4 shows three CRT displays; they are a portion of the dynamic check. 

(Note that symbols on CRT display in figure 4 are not the standard symbols defined in the 
Symbols section.) Permanent pictures of the CRT displays were obtained directly from 
the hard copy unit in the facility or from postprocessing of the plotting routine in the com- 
puter program. The plotting routine generated figure 4 by plotting the CRT display and 
adding the additional labeling on the right. 

The information output consisted mainly of calculated data preselected by the analyst 
and routed to the high-speed printer. The information could be printed for any iteration 
by activating a logical control switch. The printer is located in the proximity of the pro- 
gram control station and easily accessible to the computer operator. The output con- 
sisted of the following information: 
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APPENDIX D - Concluded 


(1) run and iteration numbers 

(2) covariance matrix of the measurement noise ^R°(N)j, its determinant, and its 

inverse 

(3) variables of state and parameters active in algorithm 

(4) nominal parameter values a? (i = 1, 2, . . p’) 

(5) calculated changes in the nominal values Aai (i = 1, 2, . . p’) 

(6) covariance matrix for the parameters in a modified form more readable to the 

analyst 


■'O!, 


Vi 


« 1«2 


'ar 


"l“p’ 


"2 V 




p’""l 




p '^2 


“p’ 


in that the standard derivations and correlation coefficients are expressed 
explicitly. 


The integration scheme that was used for the parameter estimation procedure was 
second-order Adams -Bashforth, a 1-pass integration scheme. The real-time system 
provides the option of four integration schemes: (1) second-order Runge-Kutta (2 pass), 
(2) fourth -order Runge-Kutta (4 pass), (3) second-order Adams -Moulton (2 pass), and 
(4) fourth-order Adams -Moulton (2 pass). The Adams -Bashforth scheme was obtained 
by program logic limiting scheme (3) to a 1-pass operation; this thus reduced the com- 
putation time for the integration of the equations of motion and sensitivity equations. The 
dynamic check was run by using the Adams -Bashforth scheme and scheme (2); the indica- 
tion was that the Adams -Bashforth scheme was adequate for the parameter estimation 
procedure. 
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APPENDIX E 


TEST CASES 

The test cases were for the longitudinal motion of the aircraft and included the 
effects of measurement noise on the pseudo data. 

The equations of motion are 


Cl = -qw - g sin 0 + 1 ^ V^S (Cx,o) 

(El) 

w= qu + g cos 

0 + 5 ^ 1- Czj^Be) 

(E2) 

4=,1^v2s5I 


(E3) 

0 = q 


(E4) 


where 

jD.l sin 2.5t (0 g t S tt/1.25) 

0 (tg7r/1.25) 

V = 'Ju^ + w^ 

Q!o = tan~^ ^ 

d, u 

The test cases were for different noise levels of 1, 2, 5, and 10 percent on the 
variables u, w, q, and 6 . Table I shows the known and calculated standard deviations 
of the noise for each percent level. The calculated standard deviations agreed closely 
with the known input, with an error of less than 1 percent. Table II shows the true and 
calculated parameter values and their standard deviations at each noise level. The cal- 
culated parameter values indicate convergence to within one standard deviation of the true 
values based on a fixed number of iterations. Figure 5 shows the CRT display of the con- 
verged solution and the pseudo flight data. (Note that symbols on CRT display in figure 5 
are not the standard symbols defined in the Symbols section.) 
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TABLE I.- MEASUREMENT NOISE STATISTICS USING 
MAXIMUM LIKELIHOOD ESTIMATION 


Percent 

noise 

level 

Standard deviation of - 

u 

w 

q 

e 

Known 

Calculated 

Known 

Calculated 

Known 

Calculated 

Known 

Calculated 

1 

0.05 

0.049511 

0.3 

0.29877 

0.002 

0.0019988 

0.002 

0.0019973 

2 

.10 

.099036 

.6 

.59736 

.004 

.0039965 

.004 

.0039946 

5 

.25 

.24763 

1.5 

1.4930 

.010 

.0099891 

.010 

.0099855 

10 

.50 

.49529 

3.0 

2.9854 

.020 

.019977 

.020 

.019971 
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TABLE n,- PARAMETER VALUES AND STANDARD DEVIATIONS USING MAXIMUM LIKELIHOOD ESTIMATION 


Parameter 

True 

value 

1- Percent 
noise level 

2-Percent 
noise level 

5-Percent 
noise level 

10-Percent 
noise level 

Calculated 

value 

Standard 

deviation 

Calculated 

value 

Standard 

deviation 

Calculated 

value 

standard 

deviation 

Calculated 

value 

Standard 

deviation 

^X,o 

0.112 

0.11181 

0.0001748 

0.11162 

0.0003494 

0.11105 

0.0008720 

0.11010 

0.001739 

^Z,o 

-1.29 

-1.2899 

.0007873 

-1.2898 

.001574 

-1.2895 

.003925 

-1.2892 

.007820 


-4.59 

-4.5763 

.02246 

-4.5616 

.04488 

-4.5159 

.1120 

-4.4353 

.2231 


-4.93 

-4.9332 

.04083 

-4.9331 

.08167 

-4.9313 

.2043 

-4.9286 

.4090 

Cm,o 

.0199 

.019855 

.00008130 

.019830 

.0001624 

.019747 

.0004048 

.019600 

,0008054 


-.836 

-.83457 

.001664 

-.83316 

.003326 

-.82889 

.008300 

-.82162 

.01655 

Cmq 

-32.0 

-32.102 

.1222 

-32.176 

.2446 

-32.401 

.6130 

-32.798 

1.231 

c-a. 

-3.1 

-3.1018 

.005563 

-3.1001 

.01113 

-3.0948 

.02785 

-3.0863 

.05584 


> 1 ^ 

-3 














Figure 2.- Operational diagram of Langley real-time simulation system 

for parameter estimation. 
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(a) Longitudinal motion. 

Figure 4.- Graphic display of pseudo and calculated variables 
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Figure 4.- Continued. 
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Figure 5.- Graphic display of pseudo flight and calculated 
(converged solution) longitudinal motion. 


NASA-Langley, 1972 19 L-8178 


55 


WIM/SEC) 





NATIONAL AERONAUTICS AND SPACE ADMISTRATION 
WASHINGTON. D.C. 20546 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
■ SPACE ADMINISTRATION 


OFFICIAL BUSINESS 

PENALTY FOR PRIVATE USE $300 


FIRST CLASS MAIL 



003 001 C1 U 19 720421 S00903DS 
DEPT OF THE AIR FORCE 
AF WEAPONS LAB (AFSC) 

TECH LIBRAHY/WLOL/ 

ATTN: E LOU BOWMAN, CHIEF 
KIRTLAND AFB NM 87117 


POSTMASTER: 


If Undeliverable (Section 138 
Postal Manual ) Do Not Return 


"The aeronautical and space activities of the United States shall be 
conducted so as to contribute . . . to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof.” 

— National Aeronautics and Space Act of 1958 


NASA SaENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

" CONTRACTOR REPORTS: Scientific and 

V (y^ technical information generated under a NASA 
' , ' contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs^ 
Technology Utilization Reports and 
Technology Surveys. 


Qetails on the availability of these publications may be obtained from: 

’ SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Waihinglon, D.C. 20546 



